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FOREWORD 


This  report  documents  a  FORTRAN  routine  LINOPT  for  solving  linear 
programming  problems.  Upper  and  lower  bounds  on  all  variables  are  permitted, 
and  the  dual  problem  includes  as  a  special  case  linearly-constrained  minimum 
^j-norm  problems.  Basic  theory,  the  algorithm  used,  input-output  procedures 
and  examples  of  use  are  included. 

<5 7K 

IRA  M.  BLATSTEIN 
By  direction 


r .  1  •••  •  -1  \  'n/ 


1 

i 


'.lily  Cod? 
,■  ...  r.iu'i/OF 


L 


1/2 


NSWC  TR  80-413 


CONTENTS 

Chapter  Page 

1  INTRODUCTION .  5 

2  PROBLEM  FORMAT  AND  PROGRAM  USE .  7 

3  SUPPORT  FUNCTIONS  AND  DUALITY  IN  LINEAR  PROGRAMMING .  13 

4  DUAL  SIMPLEX  METHOD .  37 

5  EXAMPLES .  21 

6  LISTING .  27 

ILLUSTRATIONS 

Figure  Page 

1  ILLUSTRATION  OF  DUAL  OPTIMAL  SOLUTIONS. 

(113*  ■  0)  . . . . .  14 

2  TEST  PROGRAM .  26 

TABLES 

Table  Page 

1  FILLING  IN  MISSING  BOUNDS  .  9 

2  CORRESPONDING  VARIABLES  .  10 


3/4 


NSWC  TR  80-413 


CHAPTER  1 
INTRODUCTION 


REPORT  ORGANIZATION 


This  report  documents  a  FORTRAN  subroutine  (with  associated  subroutines) 
called  LINOPT  for  solving  linear  programming  (LP)  problems.  It  is  divided 
into  several  chapters.  Following  the  INTRODUCTION  is  a  chapter  (PROBLEM 
FORMAT  AND  PROGRAM  USE)  explaining  the  types  of  problems  which  can  be  solved 
by  LINOPT,  some  manipulations  on  them  and  correspondences  with  program 
notation.  A  chapter  (SUPPORT  FUNCTIONS  AND  DUALITY  IN  LINEAR  PROGRAMMING) 
discusses  some  of  the  duality  concepts  behind  the  formulation  and  solution  of 
LP  problems.  The  next  chapter,  DUAL  SIMPLEX  METHOD  gives  some  information 
about  the  algorithm  used  in  the  program.  EXAMPLES  and  a  LISTING  follow. 

PROBLEM  FORMULATION 


A  rather  abstract  formulation  of  an  LP  problem  in  the  following:  Let  X  be 
a  real  vector  space  paired  with  M  under  the  bilinear  form  (inner  product) 

(p,x)  - ►  u-x,  utM,  xeX,  and  let  Y^,  i  *  1, . m  be  similarly 

paired  with  i  *  I,...,  m.  Given  ueM,  closed  convex  sets  in  Yj.,  and 
linear  transformations  A±,  :  X-*Y±f  1  •*  l,...,m, 

maximize  y-x 

subject  to  Ajx  e  Ci,  i«l,...,m. 

This  problem  in  convex  programming  has  the  following  as  its  dual  problem: 

minimized''  ^C^A^) 

*-'i*l 

Em  * 

A.|  X .  *  y 
i-1 

where  CTCi  is  the  support  function  of  the  set  Cj  —  see  Chapter  3  for 
more  on  support  functions.  For  explanatory  purposes  it  is  sufficient  to  take 
each  A^  equal  to  the  identity,  so  that  Yi  ■  X,  Aj  *  M,  i  ■  L,...,m.  Chapter  3 
treats  this  simplified  version. 


The  program  LINOPT  is  set  up  to  handle  constraints  of  the  form  u^.x  e  C^ 
where  ^  e  M  and  C^  is  a  nonempty  closed  Interval.  If  some  C^  is  a 
bounded  interval  of  nonzero  length,  the  dual  problem  has  a  nonlinear 
objective;  it  is,  however,  convex  and  piecewise  linear.  The  details  are  given 
in  Chapter  2,  where  it  is  shown  that  the  general  form  of  the  dual  objective  is 
the  sum  of  two  terms,  one  linear,  and  one  a  weighted  norm. 
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The  algorithm  used  in  the  program  is  a  form  of  the  revised  dual  simplex 
algorithm  modified  to  handle  upper  and  lower  bound  constraints.  The  inverse 
matrix  used  is  a  row-basis  inverse.  Accordingly,  the  algorithm  is  more 
efficient  on  problems  with  many  constraints.  (On  problems  with  fewer 
dependent  variables  than  independent  variables,  a  column-basis  inverse  would 
be  smaller.)  Use  of  a  row  basis  has  definite  advantages  in  modifying  a 
problem  and  then  reoptimizing. 

No  new  theory  is  involved  in  this  program.  The  dual  simplex  algorithm  can 
be  found  in  standard  linear  programming  texts. Insisting  that  all 
variables  in  the  primal  problem  have  both  upper  and  lower  bounds  makes  it 
trivial  to  find  a  dual-feasible  point  to  start  the  algorithm,  and  because  of 
the  resulting  asymmetry  between  primal  and  dual  problems,  allows  us  to  handle 
directly  (via  the  dual  problem)  certain  piecewise  linear  convex  minimization 
problems . 

The  results  from  convex  analysis  used  in  Chapter  3  can  be  found  in  greater 
generalization  and  detail  in  Rockafellar' s  book. 3 


^Hadley,  G. ,  Linear  Programming,  Addison-Wes ley ,  Reading,  1962. 

^Sinmonard,  M. ,  Linear  Programming,  Prentice-Hall,  Englewood  cliffs,  1966. 
^Rockafellar ,  R.  T. ,  Convex  Analysis,  Princeton  University  Press,  Princeton 


1970. 
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CHAPTER  2 

PROBLEM  FORMAT  AND  PROGRAM  USE 


INTRODUCTION 


LINOPT  is  programmed  to  solve  a  problem  maximizing  (or  minimizing)  a  linear 
function  subject  to  upper  and  lower  bounds  on  linear  constraint  functions. 

(The  bounds  are  equal  for  an  equality  constraint.)  The  dual  to  this  problem 
has  a  piecewise  linear  objective  function  and  linear  equality  constraints. 
Missing  bounds  in  the  primal  problem  correspond  to  sign  constraints  on  the 
dual  variables.  Such  missing  bounds  can  be  handled  by  introducing  a  penalty 
function  for  the  sign  constraint  violations  in  the  dual  problem.  An  even 
simpler  and  more  direct  interpretation  is  that  the  missing  bounds  can  be 
replaced  by  bounds  so  large  in  magnitude  that  they  are  effectively  infinite. 

PRIMAL  PROBLEM 


Maximize  x^  subject  to  <_  x^  £  k*l,...,  m+n,  where  kg  e 

(1 , . . . . ,  m+n}  Sid 

m 

Xn+i  "52  aij  Xj*  i*1 . m* 

The  objective  variable  xkg  can  be  expressed  by: 

x*o  =22.=1  cj  xj’ 

where 

c  ,  5k0j  if  kQ  1  n, 

ak0~n>j  if  kg  >  n. 

The  same  data  also  define  the  dual  problem. 

DUAL  PROBLEM 

Minimize^  ^  max  (ukbk,  ukbk)  subject  to 

Uj  Utl+i  a«  *  V  . ”• 

The  dual  objective  has  another  form  which  is  more  likely  to  be  recognized  in 
an  application: 

Zm+n  _  -  "^rmfn  /  \  __^m+n  / 

max{ukbk,  ukbv)  -  >  ,  I  \  +  5k  1  uk  +  /  .  f  bk  "  -k  »  lui 

k-1  ^k-l  1-^-2 - j  ^k-1  l - 5 - 
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Thus  the  dual  objective  contains  a  linear  term  and  a  weighted  2,^-norm  term. 

The  dual  variable  ujj  can  be  thought  of  as  a  Lagrange  multiplier  for  Che 
constraint  x^elb^,  b^]. 

MISSING  BOUNDS;  SIGN  CONSTRAINTS 


The  following  table  (TABLE  1)  shows  how  to  prepare  primal  problems  with 
missing  constraints  or  dual  problems  with  sign  constraints.  M  is  a  very  large 
positive  number.  (The  default  value  supplied  by  the  program  is  10^-®®.) 

The  first  line  of  the  table  gives  the  standard  two-sided  constraint  assumed 
by  the  program.  The  other  lines  give  the  modifications  for  unilateral  and  no 
constraints.  In  the  modified  problem  uj<  is  always  unconstrained  in  sign.  The 
dual  objective  picks  up  the  original  linear  term  when  the  sign  constraint  is 
satisfied,  and  a  penalty  term  when  the  sign  constraint  is  violated. 

When  the  program  gives  an  optimal  solution  in  which  =  +  M  for  some  k, 
the  original  problem  has  an  unbounded  solution.  If  it  has  a  finite  solution, 
the  program  will  yield  it,  and  it  will  not  depend  on  M  (unless  M  has  been  set 
so  small  that  it  interferes  with  the  "legitimate"  constraints).  The 
calculations  are  arranged  so  that  roundoff  errors  due  to  the  disparity  in 
magnitude  between  M  and  the  original  data  do  not  propagate  from  iteration  to 
iteration,  and  appear  within  an  iteration  only  if  some  x^  =  +  M. 

The  objective  variable  x^q  is  also  formally  a  constrained  variable, 
although  generally  the  constraint  will  be  -M  £xfc(£M,  i.e.  essentially  no 
constraint  at  all.  Tighter  bounds  may  at  times  be  useful.  The  constraint 
bfcg.Sxk-iM  can  be  used  to  answer  the  question:  Is  max  x^^h^?  If 
the  answer  to  this  question  is  negative,  the  constraints  are  inconsistent.  As 
soon  as  the  inconsistency  is  detected,  the  program  returns  to  the  calling 
program  without  going  on  to  calculate  the  solution  completely. 
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PROGRAM  NOTATION 


Correspondence  between  the  notation  herein  and  the  notation  of  LINOPT  is 
given  in  Table  2. 


TABLE  2  CORRESPONDING  VARIABLES 


This  TR 

LINOPT 

m 

M 

n 

N 

k0 

IOBJ 

xk 

X(K) 

uk 

U(K) 

ik 

BL(K) 

bk 

BU(K) 

aij 

A(R0W(l)  +  COL(J)) 

M  J 

BIGM 

The  FORTRAN  variable  BIGM  is  included  in  the  program  for  the  user's 
convenience.  It  supplies  a  default  value  (which  can  be  changed)  for  filling  in 
the  missing  bounds.  (The  user  must  fill  in  all  bounds,  since  there  is  no 
provision  for  keeping  track  of  missing  bounds  otherwise.) 

CONSTRAINT  COEFFICIENT  STORAGE 


The  constraint  coefficients  are  referenced  in  an  unusual  but  flexible  way. 
Row  and  column  pointer  arrays  ROW  and  COL  are  used  to  index  an  array  A.  The 
FORTRAN  standard  for  array  storage  is  by  columns  (a^,  a2i,  831,..., 
a12>  a22»  a32>*-*>  etc.)  Suppose  that  we  have  a  matrix  A  stored  in  an 
array  dimensioned  10  x  20  and  we  wish  to  study  a  problem  whose  constraint 
coefficients  form  a  submatrix  of  A,  as  in 


Vs 

V8 


We  can  set  X3  =  z-j ,  X2 


a57 

a87 


z3  •  x3 


Zt), 


a53  a59 

z7 

a83  a89 

z3 

- 

z9 

x4  =  ys>  x3  =  ys. 

ROW ( 1 )  =  5,  R0W(2)  =  8, 

C0L(1)  =*  60,  C0L( 2 )  =  20,  C0L(3)  =  80. 


COL(J)  is  set  to  the  number  of  elements  in  the  array  preceeding  the  coefficient 
column  for  X(J) —  6  x  10  for  the  7th  column,  2  x  10  for  the  3rd  and  8  x  10  for 
the  9th.  ROW(I)  then  picks  off  the  appropriate  entry  in  the  column.  We  can 
even  introduce  (by  using  the  L0CF  funct ion) columns  extraneous  to  the  array 
storing  A.  Consider  the  modified  example: 
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in  which  X3  »  W4.  COL(3)  is  set  by 

COL(3)  =  LOCF  (Bl,4) )  -  LOCF  (A(l,l)). 

The  columns  of  B  must  be  structured  the  same  way  as  the  columns  of  A  for  this 
procedure  to  work.  Instead  of  having  a  coefficient  matrix  stored  in  the  array 
in  the  usual  way  we  could  have  it  transposed,  perhaps  as  a  result  of  starting 
with  a  dual  problem.  This  corresponds  to  storage  by  rows  (not  the  FORTRAN 
standard).  A  little  reflection  indicates  that  defining  ROW  the  way  COL  is 
defined  above,  and  COL  the  way  ROW  is  defined  handles  this  storage  arrangement. 
Further  examples  are  given  in  the  program  comments. 

INPUT  AND  OUTPUT 

Input  and  output  variables  are  clearly  indicated  in  the  program  comments. 
(The  program  is  listed  in  the  section  titled  "LISTINGS".)  Arrays  are  passed  as 
formal  parameters.  Scalers  are  passed  by  using  a  labelled  common  block  /XXXLP/, 
which  must  accordingly  be  a  common  block  in  the  calling  program. 

ROUNDOFF  CONTROL 

In  the  program  there  are  three  input  variables  which  can  be  used  to  control 
roundoff  error  accummulations.  EPS  is  a  tolerance  used  in  checking  constraint 
violations.  H  is  also  used  to  zero  out  coefficients  in  the  tableau  which  have 
small  nonzero  values  (typically  for  a  CDC  60-bit  machine,  on  the  order  of 
10-i^)  which  ought  to  vanish.  Finally  constraint  violations  less  than  EPS  are 
eliminated  from  the  optimal  solution  before  returning.  For  small  problems  EPS  = 
0  is  generally  all  right. 

The  other  roundoff  -  controlling  variables  are  invert  (a  logical  variable) 
and  ITMAX.  When  INVERT  is  TRUE  the  inverse  matrix  corresponding  to  the  current 
key  K  is  calculated.  ITMAX  is  a  limit  on  the  number  of  iterations.  When  this 
limit  is  reached  control  is  returned  to  the  calling  LINOPT  again  with  INVERT  « 
TRUE.,  one  may  control  the  building  of  roundoff  error  in  the  inverse  matrix 
(which  otherwise  is  updated  by  column  operations  every  iteration).  (For  small 
problems  this  may  not  be  necessary.) 

After  obtaining  a  solution  (or  after  any  return  from  an  initial  call  to 
LINOPT)  INVERT  can  be  set  to  FALSE  for  another  call  to  LINOPT.  Certain 
modifications  to  the  problem  data  are  permissible  at  such  a  time  -  constraints 
may  be  added,  for  example.  These  modifications  are  any  for  which  the  inverse 
matrix  would  be  unchanged,  and  include  the  following  (Note:  primary  indices: 
K(l) , . . . ,K(N) ;  secondary  K(N  +  1),...,K(N  +  M)  -  see  Chapter  4.) 

ADDING  CONSTRAINTS 

M  is  increased,  new  elements  to  ROW  are  added  to  point  to  the  new  constraint 
coefficients  (which,  if  not  already  defined  should  be  stored  appropriately),  and 
new  bounds  added  to  BL  and  BU. 
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MODIFYING  CONSTRAINTS 

Bounds  for  any  secondary  variable  can  be  changed.  Inactive  bounds  for 
primary  variables  can  be  changed.  Active  bounds  for  primary  variables  can  be 
changed,  provided  the  corresponding  solution  is  also  changed;  e.g.  if  X(K(1))  31 
BU(K(1))  and  BU(K.(1 ))  Is  changed,  X(K(1))  must  be  changed  in  the  same  way. 
Constraint  coefficients  for  dependent  variables  X(N  +  1),...,  X(N  +  M)  which  are 
also  secondary  variables  can  be  changed;  (these  may  include  the  objective 
variable)  or  such  constraints  can  be  dropped,  with  appropriate  changes  to  ROW, 

M,  BL  and  BU.  (If  the  constraint  corresponding  to  ROW(I)  is  dropped,  the 
simplest  way  to  make  these  changes  is  to  set  ROW(I)  ■  ROW(M),  BL(N  +  I)  ■  BL(N  + 
M),  BU(N  +  I)  -  BL(N  +  M),  and  then  M  -  M  -  1,  so  that  the  indexing  for  X(N  +  M) 
is  changed  to  X(N  +  I)). 

Of  course,  when  LINOPT  is  recalled  with  INVERT  ■  TRUE,  any  problem  changes 
whatsoever  are  permissible. 
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CHAPTER  3 

SUPPORT  FUNCTIONS  AND  DUALITY  IN  LINEAR  PROGRAMMING 


The  support  function  Oq  of  a  convex  set  C  in  X  is  a  convex  function  defined  on 
M  by: 


ac  -  (y)  -  suPx£C  u*x 


(l) 


The  support  function  of  the  empty  set  is  -  00  everywhere.  For  nonempty  C, 

Oq  (p)  >  —  °°  and  may  take  the  value  +  in  Rockafellar's  terminology,  it  is  a 
proper  convex  function^. 


Many  of  the  formulas  of  convex  analysis  can  be  simplified  when  they  are 
restricted  to  convex  polyhedra  and  convex  polyhedral  functions.  One  such  is 
found  in  Corollary  16.4.1  of  Rockafellar's  book5,  from  part  of  which  we  can 
derive  the  following:  Let  Cl,...,  Cm  be  closed  convex  polyhedra  with  nonempty 
intersection  C  (also  a  convex  polyhedron).  Then 


mlnEiIi  °Ci  <Wi> 
subject  to  E"  »i  *  V 


(2) 


(The  general  version  of  the  corollary  is  required  if  C  is  empty.)  Rockafellar 

terms  the  operation  in  (2)  "infimal  convolution",  since  for  m  =*  2, 

c*C1nC2^  *  in£\  (aCia)  +  aQ2(MrA))  a  form  reminiscent  of  integral  convolutions. 

Formulas  (1)  and  (2)  express  Oq  (u)  as  the  common  optimal  value  of  a  pair  of  dual 
convex  programming  problems: 

Primal  Problem: 


Maximize  u*x 
subject  to  x  e  C^,  i*l,. 
Dual  problem: 

Minimize \  '  m  aci(Ui) 

A— f  i-1 

V  '  m 

subject  to/  ■  y. 

^“*i-l 


j  01  • 


4.  Ibid,  p.  24. 

5.  Ibid,  p.  146. 
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Three  solution  cases  arise: 

1.  °C  (u)  *  -  °°:  C  is  empty.  The  dual  objective  is  unbounded  below  on 

the  dual  constraint  set. 

2.  -  oo  <  CJc  (u)  <  +  °°:  C  is  not  empty.  The  optimal  value,  viz  cq  (u) ,  is 

attained  in  both  problems. 

3.  0C  (y)  >  +  oo;  C  is  not  empty.  The  primal  objective  is  unbounded 

above  on  the  primal  constraint  set.  The  dual  objective 
is  +  oo  everywhere  on  the  dual  constraint  set,  (i.e. 
aCi  (ui)  <  00 ,  i  *  l,...,m  =>  EjSi  Ui  i  u. 

Suppose  that  °C  (u)  is  finite,  and  let  x*  solve  the  primal,  Ui*,  i*l,...,m,  the 
dual.  Then 

aC  (u)  *  U’X*  (primal  optimality) 

3  Ui*)‘x*  (dual  constraint) 

"  Ei-1  <^i*  '  x*> 

<  (Uj*)  (primal  constraints) 

3  aC  (u)  (dual  optimality) 

It  follows  that 

Ui*  ’  x*  3  %  (Ui*),  i“l,  •  •  .  ,m,  (3) 

or  that  Ui*  supports  Ci  at  x*:  Ui*  is  an  outer  normal  to  Ci  at  x*.  The  formula 
is  one  way  of  expressing  complementary  slackness,  since  if  x*  Z  int  Ci  then  Ui*  =  0, 
while 


Figure  1  ILLUSTRATION  OF  DUAL  OPTIMAL  SOLUTIONS.  (U3*  3  0) 
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if  Ui*  ■  0,  then  Che  consCrainC  x  e  Cf  can  be  dropped  withouc  alcering  the 
soluCion.  The  set  of  constraints  for  which  yj*  4  0  are  active  (binding) at  the 
optimum.  (See  Figure  1).  There  may  be  multiple  solutions,  each  with  a 
different  set  of  active  constraints. 

The  usefulness  of  (2)  hinges  on  picking  the  constraint  sets  C^,  i=l,...,m,  to 
be  simple  enough  to  permit  easy  evaluation  of  their  support  functions.  Any  convex 
polyhedron  can  be  expressed  as  the  intersection  of  half -spaces,  and  any  half -space 
can  be  defined  by  a  linear  inequality.  Let  H:  =>  (x:  y*x  £  b},  where  y  4-  0.  Then 


crH  (v) 


!ub  if  v 
+  00  othi 


uy  and  u  >  0, 


otherwise. 


In  terms  entirely  of  hyperplane  constraints  the  primal  and  dual  problems  become 
Primal  problem: 

Maximize  y*x 

subject  to  y^  -  x  _<  b^,  i*l,...,m 


Dual  problem: 


Minimize 


\  m 

Lj 


subject  toV'  m  ut  y^_  ■  y, 

1 

u-l  >  0,  i«l, 


The  sign  constraints  on  u^  avoid  infinite  values  of  the  dual  objective  and  keep 
it  linear  so  that  both  problems  consist  of  optimizing  linear  functionals  subject 
to  linear  constraints.  Alternatively  we  could  omit  the  sign  constraints  and  keep 
the  formulation  in  terms  of  support  functions.  This  pair  of  problems  also 
illustrates  the  more  general  duality  relationship  cited  in  the  INTRODUCTION.  Let 
A^:  X-*R  be  defined  by  Aj_X:  *  y^  «  X.  Then  Aj*:  R-*M  and  A^*u  »  uy^.  Moreover, 


<,<-.b1](u>  ■ 


ub^  if  u  >  0, 
+  oo  if  u  <  0. 


Thus  the  problems  are  expressible  as: 

Primal  problem: 

Maximize  y*x 

subject  to  A^  x  e  (-“>,  bf] ,  i*l,...,m 


Dual  problem: 


Mlnimize^r^  m 

Wi-1 


(-“.bil 


subject  to  e,:, 


A.j*  Uj_  *  y 
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We  may  have  two-sided  constraints  such  as 

b  <  yx  <_  b 

(where  b  <_  b)  defining  a  set  S  in  X,  which  can  be  replaced  by  the  pair  of 
constraints 


P'x  £  b, 
-y  x  <  -b. 


since 

!max  {ub,  ub}  if  v  *  up, 

+  00  otherwise, 

and  max  (ub,  ub)  is  nonlinear  in  u  (unless  b  ■  b) ,  the  dual  of  a  problem  with 
two-sided  constraints  is  nonlinear.  Of  course  it  is  easy  to  relate  the  two-sided 
and  one-sided  versions  by  using  (2).  Thus  if  u  is  the  dual  variable  for  the 
constraint  b  <_  yx  <_  b,  u+  for  yx  £  b  and  u~  for  -yx  <_  -b, 

then  u  *  u+  -  u"  (4) 

while  if  b  ^  b  either  u+  or  u“  vani*shes  (at  the  solution),  so  that 

u+  51  max  {u,0},  -  u~  ■  min  {u,0}.  (5) 

When  b  -  b,  only  the  difference  u  *  u+  -  u~  is  determined.  An  alternate 
viewpoint  in  this_case  is  that  u  is  the  dual  variable  for  the  linear  equality 
constraint  yx  *  b. 

The  pivoting  operations  of  the  dual  simplex  method  can  be  thought  of  as 
substituting  one  hyperplane  bounding  a  half -space  for  another,  and  consequently 
are  better  suited  for  the  formulation  in  terms  of  one-sided  constraints.  The 
relations  (4)  and  (5)  and  some  sign  bookkeeping  then  make  it  easy  to  apply  the 
method  to  two-sided  constraints.  Explanations  of  the  method  without  the  sign 
manipulations  are  more  transparent.  Accordingly  in  the  next  section  only  one-sided 
constraints  are  considered. 
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CHAPTER  4 

DUAL  SIMPLEX  METHOD 

In  the  previous  sect '  \n  no  particular  coordinates  were  used  on  X.  Most  LP 
problems  encountered  are  expressed  in  terms  of  coordinates  with  respect  to  some 
particular  basis  for  X,  the  coordinates  then  forming  a  set  of  independent 
variables,  (as  indicated  previously,  LINOPT  assumes  such  a  formulation.)  Thus, 
with  one-sided  constraints,  we  get  a  pair  of  problems  like  the  following,  in 
which  we  assume  that  the  n  independent  variables  xj ,  jeJ  are  included  in  the  nri-n 
constrained  variables  x^,  keK;  i.e.  JCK. 

Primal  problem: 

Maximize  Xq  :  *  eJ  cj  xj 

subject  to  x^:  »  £jej  a^j  Xj  £  b^,  keK. 

(Note  that  a^  *  5  for  keJ.) 

Dual  problem: 

Minimize  E  keK  uk  bk 
subj  ect  to  I  k£K  u^  a^j  *  Cj  ,  j  eJ 
uk  >_  0,  keK. 

Note  that  the  n  equations  relating  the  dual  variables  u^,  keK  can  be  written 
explicitly  for  Uj : 

2  ieK~J  ui  aij  +  uj  "  cj 

Thus  in  the  dual  problem,  u^,  ieK~J,  are  independent  and  Uj ,  jej,  are  dependent. 
Given  some  other  subset  J’  of  K  for  which  Xj ,  jeJ’,  are  linearly  independent, 
we  can  transform  the  constraint  relations  so  that  xj ,  jeJ'  are  the  independent 
variables  through  which  the  primal  problem  is  phrased.  Given  such  as  index  set 
J'  we  can  define  a  corresponding  basic  solution .  For  J  the  definition  of  a 
basic  solution  is  obtained  by  setting  the  independent  variables  Xj ,  j£j,  and 
Uj,  icKM,  to  their  bounds  and  satisfying  the  constraint  relations  among  the 
variables;  x^  and  u^  are  the  values  of  x^,  at  the  basic  solution. 

xj  *  bj ,  j  £ J  (primary  primal  variables) 

xi  *  Z  jej  a.,  b* ,  i£K~J  (secondary  primal  variables) 

j  .  j,,  jeJ  ^  (secondary  dual  variables) 

’  (primary  dual  variables) 

Ui  -  0,  i£K~J 
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The  terms  primary  and  secondary  have  been  introduced  instead  of  independent  and 
dependent  because  one  may  wish  to  refer  to  the  primary  variables  of  the  problem 
as  initially  formulated  as  the  independent  variables.  The  split  between  primary 
and  secondary  depends  on  the  set  J  and  changes  with  it.  The  secondary  primal 
indices  are  usually  called  basic  indices  in  linear  programming  texts,  since  they 
correspond  to  the  indices  for  a  column  basis  for  the  constraint  matrix.  This 
terminology  is  a  little  inappropriate  here,  since  LINOPT  makes  use  of  a  row 
basis  corresponding  to  the  complementary  set  of  indices  -  the  dual  basic  indices 
in  the  usual  description.  The  use  of  "basic"  in  this  sense  is  avoided  here  to 
prevent  confusion. 

Furthermore,  in  a  problem  with  two-sided  constraints  the  basic  indices  refer 
to  the  indexing  of  the  equivalent  one-sided  problem,  not  the  indexing  of  the 
two-sided  problems,  so  that  the  basic  variables  for  J  would  be  x^  and  -x^ 
for  keK-^J  and  either  xj  or  -x;  (but  not  both)  for  jej.  Alternatively  we  may 
retain  the  "primary /secondary"  notation  and  supplement  it  with  some  way  of 
indicating  whether  a  primary  variable  is  at  its  upper  or  its  lower  bound.  (The 
program  simply  checks  the  solution  value  against  the  bounding  values. 

In  a  basic  solution  the  primary  variables  satisfy  the  constraints  placed  on 
them  if  all  primal  variables  satisfy  the  constraints,  the  basic  solution  and  J 
are  primal-feasible.  If  the  dual  constraints  are  satisfied,  the  basic  solution 
is  dual-feasible.  A  basic  solution  which  is  both  primal-  and  dual-feasible  is 
optimal.  At  a  basic  solution  both  primal  and  dual  objective  variables  have  the 
value  ^ j£ j  Cjbj.  The  criterion  for  dual  feasibility  is  simply  that 
cj  >0,  jeJ. 

The  transformation  of  the  constraint  coefficients  accompanying  a  change  from 
one  set  of  primary  variables  to  another  can  be  performed  explicitly  when  needed, 
or  it  can  be  expressed  in  terms  of  a  nonsingular  matrix  relating  the  variables. 

There  are  two  ways  of  doing  this.  Let  xj,  jeJ  and  xj*,  j'  eJ*  be  two 
sets  of  primary  primal  variables.  Set  I  *  K~J,  I'  =K~J' ,  and  let  xj  be  a 
column  vector  whose  entries  are  xj,  j ej,  etc.  Using  matrix  notation  the  two 
ways  can  be  described  as  follows; 

l.  Solution  for  xj < 

*  Axj  (A  is  ra  x  n) 

Rearrange  columns  to  give: 

B  xj '  *  R  xj> 

with  B  a  nonsingular  m  x  m  submatrix  of  [I  -A].  Then  xj  •  =  B~^-R  xj' . 

The  columns  making  up  B  are  a  basis  for  the  space  spanned  by  the  columns  of 
[I  -A] 

Applied  to  the  dual: 

UXA  +  Uj  =  C 
UT'S  +  Uj'D  -  C  _ 

Uj '  SD“1  +  Uj<  *  C,  where  C  *  C 
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The  rows  of  D  form  a  basis  for  Che  space  spanned  by  Che  rows  of 

2.  Solution  for  xj  and  subsCiCuCion: 


Rearrange  rows  Co  give: 


(S  and  D  are  Che  same  as  above.  B~^R  *  SD“^.) 


For  che  dual: 


IUX  Uj]  =  Ux I I  -A]  +  t0  C] 

[UX «  Uj*]  »  UX[B  -R]  +  tcx •  Cj'l 

[ux<  Uj>]  -Ux'tl  -B"lR]  +  (0  ?] 

where  S’  »  Cj*  +  CX'  B“^R  and  is  Che  same  as  before.  The  inverse  matrix 
is  the  produce  of  elementary  row  operations;  D“1  is  Che  produce  of 
elementary  column  operations.  Either  one  may  be  used  to  keep  track  of 
changes.  LINOPT  uses  D-*-  and  generates  coefficients  and  solutions  from  the 
original  constraint  coefficients  by: 


When  a  basic  solution  is  changed,  D-^  is  updated  by  column  operations.  (We 
have  ignored  xq :  assumme  that  OeK  so  that  c  is  a  row  of[V|\ 

LlJ/ 

The  dual  simplex  algorithm  works  with  dual-feasible  basic  solutions.  Given  a 
set  J  defining  primary  primal  variables,  J  is  altered  by  adding  an  index  not  in 
J  and  dropping  an  index  in  J:  a  secondary  variable  replaces  a  primary 
variable.  The  resulting  changes  in  the  constraint  coefficients  can  be 
accomplished  by  Gauss-jordan  pivoting.  The  indices  entering  and  leaving  J  are 
chosen  in  such  a  way  that  dual-feasibility  is  maintained  and  the  objective 
function  value  does  not  increase.  Proof  of  convergence  can  be  found  in  any 
linear  programming  text.^  The  procedure  is: 

1.  Identify  violated  constraints:  I'  *  {i£K~J:  XX  >  bj)  . 

2.  If  I'  is  empty,  stop:  basic  solution  is  optimal. 

3.  Pick  (by  some  heuristic)  tel'. 

4.  Identify  constraints  which  can  be  dropped  without  being  violated  when  x^ 

is  set  to  b^:  J' :  =  { jeJ:  a^j  >o}  ? 

5.  If  J'  is  empty,  stop:  constraints  are  inconsistent. 


For  example,  Hadley,  op.  cit, 


or  Simmonard,  op.  cit. 
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6.  Identify  subset  of  J'  corresponding  to  constraints  for  which  dual  feasibility  is 
maintained  when  dropped:  J"  ■  {jeJ1 


c3/akj 


^min  (cj/a^j)}. 

7.  Pick  (while  applying  anticycling  criterion,  if  desired)  £eJ". 

8.  Update  solution:  J:  *  JU(k}  ~  {£}.  Update  inverse  D-l  by  column 
operations.  Calculate  new  basic  solution. 

9.  Go  to  1. 


The  coefficients  a^j  and  Cj  are  those  corresponding  to  the  current  index  set 
J,  not  the  original  one  in  terms  of  which  the  problem  is  phrased,  and  are 
calculated  by  post-multiplying  an  original  constraint  matrix  row  (or  objective 
row)  by  D“l. 

In  the  program  the  heuristic  used  in  step  3  is  to  pick  the  most  violated 
constraint.  In  step  7  a  tie  for  4  is  broken  randomly,  a  procedure  which  prevents 
cycling  almost  surely. 
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CHAPTER  5 
EXAMPLES 


TEST  PROGRAM 

A  small  test  program  to  run  the  following  examples  is  listed  in  Figure  2. 

The  lines  between  the  call  to  LINOPT  and  the  call  to  TABLO  merely  do  some 
cosmetic  surgery  on  the  output,  replacing  quantities  near  M  in  magnitude 
(actually  those  >/M)  by  +  R  (machine  infinities).  Three  examples  are  given, 
with  NAMELIST  inputs  and  the  outputs  from  the  program.  Note  that  in  all  three, 
ROW  and  COL  are  defined  to  correspond  to  storage  by  rows,  and  the  columns  of  the 
array  A  contain  the  rows  of  the  constraint  coefficient  matrix. 

EXAMPLE  1 

This  example  is  essentially  the  problem  discussed  in  Section  1-3  of 
Hadley^,  with  slightly  modified  coefficients. 

Maximize  Xg;  =  5.0  xi  *  7.6  x2  +  8.0  x3  +  4.0  x^ 

subject  to  xx  >o,  x2 ^0,  x3_>0,  x4^0» 

x5:  -  1.5  x!  +  1.2  x2  +  2.4  x3  +  1.2  x4_<  2100, 

X6:  a  1.0  xi  +  4.5  x2  ♦  1.0  x3  +  3.0  X4  <  8000, 

x7:  =  1.5  x^  3.0  x2  +  3.6  x3  +  1.0  x4  <  5000. 
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NAMELIST  INPUT: 

$IN 

IOBJ  -  8 
M  ■  4, 

N  -  4, 

MIN  -  .FALSE., 

INVERT  »  .TRUE., 

ITMAX  -  1000, 

EPS  -  0., 

ROW  »  0,  10,  20,  30, 

COL  -  1,  2,  3,  4, 

A(1 ,  1)  -  1.5,  1.2,  2.4,  1.2, 

A(l,  2)  =■  1.0,  4.5,  1.0,  3.0, 

A(1 ,  3)  *  1.5,  3.0,  3.6,  1.0, 

A(1 ,  4)  *  5.0,  7.6,  8.0,  4.0, 

BL  -  4*0.,  4*-l.E100, 

BU  *  4*1.E100,  2100.,  8000.,  5000.,  1.E100, 

K  -  1,  2,  3,  4,  5,  6,  7,  8, 

X  =>  4*0., 

SEND 

(Since  the  objective  coefficients  Cj  are  all  nonzero  for  the  initial  tableau, 
it  is  not  really  necessary  to  preset  xi ,  X£,  X3  and  X4,  as  Example  2  will  show). 

Output : 


TABLEAU 


I 

BL(1) 

x(I) 

BU(I) 

T(I,  1)  T(I 

3) 

T(I,  5) 

T(I,  7) 

1 

0.000 

0.000 

R 

1.000 

0.000 

0.000 

0.000 

2 

0.000 

1625.000 

R 

-.125 

-.800 

-.417 

.500 

3 

0.000 

0.000 

R 

0.000 

1.000 

0.000 

0.000 

4 

0.000 

125.000 

R 

-1.125 

-1.200 

1.250 

-.500 

5 

-R 

2100.000 

2100.000 

0.000 

.000 

1.000 

-.000 

6 

-R 

7687.500 

8000.000 

-2.938 

-6.200 

1.875 

.750 

7 

-R 

5000.000 

5000.000 

-.000 

.000 

.000 

1.000 

8 

-R 

12850.000 

R 

-.450 

-2.880 

1.833 

1.800 

The  objective  variable. 

x(8),  is 

to  be  maximized. 

ITER 

=  10,  IERR  » 

0. 

The  tableau  gives  information  about  the  primary  and  secondary  variables  at  the 
final  iteration.  The  primary  variables  are  xj,  X3,  X5  and  X7.  The  rows  of  the 
tableau  give  the  coefficients  of  the  variables  expressed  in  terms  of  the 
primary  variables.  Thus,  X4  --1.125  x^  -1.2  X3  +  1.125  X5  -.5  xy. 

The  dual  variables  are  not  printed  explicitly  but  the  nonzero  ones  can  be 
obtained  from  the  objective  row:  uj  *  -.450,  U3  =*  -2.880,  U5  =■  1.833, 

U7  *  1.800.  (For  a  minimization  problem,  these  should  be  negated.) 
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EXAMPLE  2 

This  problem  also  cooes  from  Hadley®.  It  is  his  Problem  8-5. 

Minimize  X4:  ■  3X]_  -2x2  +  4x3 

subject  to  X1  *2  _>0,  X3>.0, 

X5:  *  3xi  +  5x2  +  4x3  >.7, 

xg:  *  6x1  +  X2  +  3x3  >_4, 

X7 :  *  7  xj  -  2X2  -  X3  <_  10 , 
xg:  *  xi  -2x2  +  5x3  >.3, 
xg:  »  4xi  +  7X2  -  2x3  >2. 

Input: 

$IN 

I OBJ  *  4, 

M  *  6, 

N  »  3, 

MIN  »  .  TRUE., 

INVERT  =  .  TRUE., 

UMAX  -  1000, 

EPS  «  0., 

ROW  =  0,  10,  20,  30,  40,  50, 

COL  -  l,  2,  3, 

A(l,  1)  *  3.,  -2.,  4., 

Ad,  2)  =■  3.,  5.,  4., 

Ad,  3)  -  6.,  1.,  3., 

Ad,  4)  *  7.,  -2.,  -1., 

A(l,  5)  -  1.,  -2.,  5., 

Ad,  6)  -  4.,  7.,  -2., 

BL  »  3*0.,  -I.E10Q,  7.,  4.,  -1.E100,  3.,  2., 

BU  *  6*l.E100,  10.,  2*1.E100, 

K  «  l, 2, 3, 4, 5, 6, 7, 8, 9, 

Send 


Ou  tpu  t : 

TABLEAU 


I 

BL(I) 

X(t) 

BU(I ) 

T(I,  l) 

T(1 ,  5) 

T(I ,  8) 

l 

0.000 

0.000 

R 

1.000 

1.000 

0.000 

2 

0.000 

R 

R 

-.333 

.152 

-.121 

3 

0.000 

R 

R 

-.333 

.061 

.152 

4 

-R 

-R 

R 

2.333 

-.061 

.848 

5 

7.000 

R 

R 

-.000 

1.000 

.000 

6 

4.000 

R 

R 

4.667 

.333 

.333 

7 

-R 

-R 

10.000 

8.000 

-.364 

.091 

8 

3.000 

3.000 

R 

.000 

0.000 

1.000 

9 

2.000 

R 

R 

2.333 

.939 

-1.152 

THE  OBJECTIVE  VARIABLE,  X(  4),  IS  TO  BE  MINIMIZED. 
ITER  -  3  I  ERR  =*  0 

8.  Ibid.,  p.  267. 
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This  problem  has  an  unbounded  soLution:  X5  =  +  00  and  the  objective  value 
is  -  00. 

EXAMPLE  3 


This  example  illustrates  the  solution  of  a  dual  problem. 

Minimize  £  juk| 

subject  to£k„3  uk  3  l» 

“^k*3  ^  uk  = 

(The  indexing  starts  at  3  for  convenience.)  Since  there  is  no  unit  matrix  in 
Che  consCrainC  coefficient  matrix,  we  add  artificial  variables  u^  and  U2 , 
which  must  vanish  at  Che  solution: 

U1  +  Ek=3u*  =  l* 
u2  +  E^k  3  l* 

Noting  that  | |  *  max  {-0^,  ,  we  can  transform  to  the  primal  problem: 

Maximize  xio:  a  X£  +  X2 

subject  to  -lix£il  ,  i  =  3 . 9,  where 


x3; 

=  xl 

+ 

3x2, 

*4  = 

3  X1 

+  4x2 

x5: 

3  xL 

5x2  , 

3  X1 

+  6x2 

x7; 
*9  • 

=  xl 
3  X1 

+ 

+ 

7X2, 

9x2. 

x8 : 

3  X1 

♦  8*2 

The  variables  and  X2 ,  dual  to  artificial  variables,  are  not  constrained 
directly. 

Input: 

$IN 

IOBJ  =  10, 

M  *  8, 

N  »  2, 

MIN  *  .FALSE., 

INVERT  ■  .TRUE., 

ITMAX  -  1000, 

EPS  »  0.  , 

RO'4  *  0,  10,  20,  30,  40,  50,  60,  70, 

COL  -  l,  2, 

Ad,  l)  *  1.,  3., 

A(l,  2)  -  l.,  4., 

A(l,  3)  -  l.,  5., 

A(l,  4)  *  l.,  6., 
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A(l,  5)  =*  l.,  7., 

A(  1 ,  6)  =■  1.,  8., 

A(l,  7)  =>  1.,  9., 

A( 1 ,  8)  »  1.,  1., 

BL  -  2*-l.ElOO,  7*-l. ,  -1.E100, 
BU  -  2*1 .ElOO ,  7*-l. ,  I.E100. 

X  -  1,2,3,4,5,6,7,8,9,10, 

X  *  2*1. ElOO, 

$END 


Ou  tpu  t : 


I 

BL(I ) 

X(I) 

BU(I ) 

T(I 

.  3) 

T(I,  9) 

l 

-R 

2.000 

R 

1.500 

-.500 

2 

-R 

-.333 

R 

-.167 

.167 

3 

-1.000 

1.000 

1.000 

1.000 

-.000 

4 

-1.000 

.667 

1.000 

.833 

.167 

5 

-1.000 

.333 

1.000 

.667 

.333 

6 

-1.000 

.000 

1.000 

.500 

.500 

7 

-1.000 

-.333 

1.000 

.333 

.667 

8 

-1.000 

-.667 

1.000 

.167 

.833 

9 

-1.000 

-1.000 

1.000 

.000 

1.000 

10 

-R 

1.667 

R 

1.333 

-.333 

The 

solution 

is  obtained 

l  from  row 

10: 

u3  =« 

1.333,  u9  =  - 

us  a0  «  The  minimal  value  is  1.667.  (Obviously,  the  exact  solution  has  u3 
*  4/3,  u9  *  -1/3.) 


25 


NSWC  TR  80-413 


cNcsnvf^cs^^xOfsoo^'-'rocNncMfn 


-3 

HWCrfWHUUUJHHHH  —  X  ^  U  U  U 

cozzzcozzzcocococoaxoazzz 

U  H  M  n •  C&3  *—  »— i  »— iUU'jJUO^OOmhh 


Nn>jMn»j'Dr>.Qoa'0'^Nf<''Ji^^r^co 
i-l  r-4  HHp-H(SNcS(NNNNNN 


HHHUUlUHHHHHHHhHHHHH 

COWCflZZZWMMMWWWWWWWWW 


a* 

1 

1 

H 

•  x- 

1 

H 

m 

1 

ca 

z 

CO 

ca 

6- 

1 

f 

*>»✓ 

w 

z 

>  r 

os 

X 

1 

* 

CJ) 

• 

» 

3 

1 

1 

• 

O  1 

CxJ 

1 

xJ 

»-4 

W 

-3 

II 

n 

O 

1 

H 

1 

2* 

• 

*  t 

to  i 

> 

• 

1 

o 

ao 

5 

O 

O 

1 

04 

1 

q 

/— N 

Z  cu 

Z  1 

z 

1 

u 

i 

HI 

oS 

H 

* 

1 

O 

1 

as 

o 

z 

i 

XX 

00 

1 

fiC 

OJ 

Oi 

1 

Z 

» 

CM 

• 

1 

2e 

z 

U3 

• 

/-N 

o 

H 

1 

xx 

1 

• 

V-^ 

-I  -  £ 

H  l 

• 

•» 

1 

O 

o 

Z  J 

cc 

XX 

o 

XX 

Z 

1 

x4 

I 

/•% 

3 

CQ  i— 

,  o 

aS  t 

z 

xJ 

1 

oS 

a.  • 

cs 

\_/ 

OS 

tt4 

1 

1 

O 

O  o, 

1  *-X 

O 3  1 

xx 

CO 

t 

•* 

CO 

z  ^ 

• 

X 

•> 

z 

1 

a! 

1 

xU 

* 

vJ 

xx  Z 

i  CO 

>  1 

z 

1 

< 

x-» 

/■v 

'X' 

to 

< 

L* 

1 

O 

1 

* 

o 

Z  1 

1 

w 

H 

•  w 

XX 

CO 

1 

c* 

1 

O 

O 

CJ 

.  ^ 

XX  1 

* 

-3 

1 

H 

OS 

-*  J 

0 

CO 

z 

<r 

1 

1 

CM 

04  CL. 

>  o* 

1 

"x* 

z 

O 

1 

04 

cr 

=> 

< 

XX 

o 

2j 

O 

r% 

co 

1 

1 

5 

1 

f 

%-/ 

< 

w 

X 

2* 

!ri 

Z  1 

z 

XX 

O 

1 

1 

O 

z 

CO 

II  W 

CO 

H 

z 

00 

cn 

w 

I 

os 

1 

q 

X  X 

:  x 

t— «  1 

■■xh 

Z 

! 

Z 

XX 

II 

h-*  bu 

Cju 

Ct4 

o 

< 

•» 

rs* 

P 

1 

o 

1 

z 

z 

os 

X  X 

:  x 

Z  1 

1 

XX 

J 

XX 

O 

H 

£ 

1 

g 

1 

o 

o 

'v.  ^ 

. 

1 

.  H 

* 

q 

1 

CM 

o 

H 

P 

H 

s 

1 

oS 

1 

hx 

XX 

as 

q  1 

CO 

OS 

1 

q 

z 

CM 

J 

Z 

z 

CO 

< 

1 

flu 

1 

CO 

CO 

w 

z  z 

;  z 

X-* 

CO 

\ 

< 

S 

Xx 

XX 

M 

OS 

1 

1 

z 

z 

o 

o  o 

i  o 

o  i 

J 

o 

* 

1 

CO 

<2 

h-4 

q 

OS 

OS 

H 

O 

1 

H 

t 

CsJ 

w 

01 

x  x 

X 

xx  l 

03 

XX 

< 

1 

os 

o 

CO 

2 

u 

o* 

04 

o 

1 

CO 

1 

s 

z 

p 

X  X 

X 

O  1 

£ 

1 

y 

QS 

1 

CO 

t 

I-* 

xx 

z 

o  o 

i  o 

O  I 

< 

1 

2 

o« 

1 

1 

H 

1 

1 

Q 

Q 

>— < 

u  u 

1  u 

xj  1 

1 

z 

XX 

CM 

1 

1 

£*J 

£  1 

►-  (i.  w 

e>  z  z 

a  h.  O 


J  ^  X 
aa  xx 

w  <N 
Z  3! 
^  3J  U 


5^  OS 

•  os 

as  r»> 

:j>  - 

CO 

C*3  S 

*  O 

X  *-< 


O  O  U  O  O  U  CJ 


t 


o  o  o  o  o  o  o  o  o  o  o  o  ra  r;  o  O  o  c*  r>  o  o  o  o  o  o  o 


NSWC  TR  80-413 


CHAPTER  6 
LISTING 


C 

C 

C 

c 

c 

c 

c 

c 

c 

{' 

c 

c 

c 


SUBROUTINE  L INQPT ( A . ROU, COL , BL ,30 ,  K  ,X ,U . E , SCR )  LI SORT 

LINE 

. — . - . LINE 

LINE 

LINEAR  PROGRAMING  by  THE  DUAL  SIMPLEX  ALGORITHM  LINOPT 

LINOPT 

LINOPT 

PROBLEM —  LINOPT 

LINOPT 

MINIMIZE  OR  MAXIMIZE  XilOBJ)  SUBJECT  TO  LINOPT 

LINOPT 

X (N+I )  =  SUMi J  =  A<RQU(I)+C0l( J) )  ■*  X'J),  LINOPT 

I  =  i , . . . , M ,  LINOPT 

LINOPT 

8L(  J)  .Lt.  X(J)  .IE.  SU( J; ,  J  =  1,...,N+rt.  LINOPT 

LINOPT 

.  LINOPT 

LINOPT 

FURTHER  DOCUMENTATION  AMD  EXAMPLES  OF  USE  CAM  BE  FOUND  IN —  LINOPT 

LINOPT 

NSUC  TS  80— A 1 3,  LINOPT,  A  FORTRAN  ROUTINE  FOR  SOLOING  LINQPT 

LINEAR  PROGRAMMING  PROBLEMS,  BY  J.U.  U INGATE.  LINOPT 

LINOPT 

- . -  LINOPT 

LINOPT 

ARRAYS  ARE  PASSED  AS  FORMAL  PARAMETERS,  SIMPLE  VARIABLES  AS  LINOPT 

ELEMENTS  OF  THE  COMMON  BLOCK  /XXXLFV.  LINOPT 

LINOPT 

INPUTS— THE  FOLLOW  I  T!G  VARIABLES  AND  ARRAYS  MUST  BE  DEFINED  ON  LINOPT 
ENTRY  LINOPT 

LINOPT 

I03J  INDEX  OF  THE  OBJECTIVE  VARIABLE,  (INTEGER)  LINOPT 

(NOTE  THAT  XCIOBJ)  IS  ALSO  CONSIDERED  AS  LINOPT 

A  CONSTRAINED  VARIABLE.)  LINOPT 


LINOPT 

M  NUMBER  OF  DEPENDENT  VARIABLES.  (INTEGER)  LINOPT 

H  NUMBER  OF  INDEPENDENT  VARIABLES.  (INTEGER)  LINOPT 

LINOPT 

r<  l  H  .TRUE.  FDR  MINIMIZATION,  (LOGICAL)  LINOPT 

•  FALSE.  FOR  iiAXIr'IZATION.  LINOPT 

LINOPT 


S3 

34 

35 

36 
3? 


-0 


27 
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C 

C 

C 

c 

c 

c 

c 

r. 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

G 

c 

c 

c 

c 

c 

c 

c 

c 


invest 


I T  ri  A  i 


.TRUE.  IF  THE  INVERSE  MATRIX  E  iL SGICAL ! 

IS  TO  3£  CALCULATED. 

•FALSE.  IF  E  IS  ALREADY  SET  TO  THE  INVERSE  FOR 
THE  FASTS  DEFINED  BY  K.  (FOR  REOPTIhlZATIGM, 

INVEST  SHOULD  BE  .FALSE.  UNLESS  REINVERSION 
IS  DESIRED.) 

MAXIMUM  NUMBER  0*  ITERATION'S  ALL'.  .ED.  ( i.'-TEG'ER) 

CONTROL  IS  RETURNED  TO  THE  CALLING 
PROGRAM  AFTER  UMAX  ITERATIONS. 

ZERO  TOLERANCE.  CONSTRAINT  VIOLATIONS  (REAL) 

OR  TABLEAU  ENTRIES  .LE.  EPS  IN  MAGNITUDE 
ARE  TREATED  AS  ZERO. 

ARRAY  CONTAINING  THE  COEFFICIENT  MATRIX.  (REAL ) 

ROW  INDEX  ARRAY.  (INTEGER) 

COLUMN  INDEX  ARRAY.  (INTEGER) 

THE  COEFFICIENT  OF  XI J)  IN  THE  EQUATION 
FOR  X(  N+I  i  IS  AIRQUdJ+COLU)). 

EITHER  (CASE  1) 

ONE  HAS  VECTORS  ARGU1 , . . . , ASGUM  UITH 
ARQUKCQL!  J) )  THE  COEFFICIENT  OF  X(J) 

IN  THE  EQUATION  FOR  XiN H>,  IN  UHICH  CASE 
ROU(I)  =  LOCF ( AROUI 5  -  LOCF(A),  I  = 

OR  (CASE  2) 

ONE  HAS  VECTORS  ACOLt . ACCLN  UITH 

ACOLJ(ROUiI))  THE  COEFFICIENT  OF  X(J) 

IN  THE  EQUATION  FOR  X(N+I>,  IN  UHICH  CASE 

COL(J)  =  LOCF(fiCCLJ)  -  LOCF! A) ,  J  =  ( . r(. 

(E.G.  IF  A  IS  DIMENSIONED  FOR  MM  ROUS 
AND  THE  COEFFICIENT  MATRIX  IS  STORED  IN 
THE  FIRST  M  ROUS  AND  N  COLUMNS  OF  A  (CASE  2>r 
ROU(I)  =1,  I  = 

COL(J)  =  ( J-l  )*MM,  J  = 

WHILE  IF  THE  COEFFICIENT  MATRIX  IS  STORED 
TRANSPOSED  IN  THE  FIRST  N  ROUS  AND  M  COLUMNS 
(CASE  !  ) , 

ROU(I)  =  (I-mrffi,  I  =  1 . a, 

COL C J)  =  J,  J  *  1,...,N. 

ROU  AND  COL  MAY  BE  PERMUTED  IN  ANY  CONVENIENT 

UAY.) 

ARRAY  OF  LOWER  BOUNDS.  (REAL) 

ARRAY  OF  UPPER  BOUNDS.  (REA,.) 

BASIC  SOLUTION  KEY.  (INTEGER) 

K,  IN  CONJUNCTION  UITH  X,  SPECIFIES  A 
PARTICULAR  BASIC  SOLUTION.  THE  EQUATIONS 

RELATING  X(Nd),  I  =  1 . M  TO 

X  (  J  / ,  •  =  t  .  •  • . , N  *  1  H c  ;  S T I  '■  1 0 U A T  i •. N 3  .* 

CAN  BE  SOLVED  FOR  VARIOUS  CCMBI.-A  '  I  j  '•  S  L." 


LlriflPT 
LlnOF  >" 
LI, SOFT 
LI  i*0?  T 
LI, '(OPT 
LINOPT 
LIHOPT 

l : no? r 

LINOPT 
LINOPT 
LINOPT 
LINOPT 
LINOPT 
LINOPT 
LINOPT 
LINOPT 
LINOPT 
LINOPT 
LINOPT 
LINOPT 
LINOPT 
LINOPT 
LINOPT 
LINOPT 
LINuPT 
LINOPT 
LIHOPT 
LINOPT 
LINOPT 
LINOPT 
LINOPT 
LINOPT 
LINOPT 
LINOPT 
LINOPT 
LINOPT 
LINOPT 
LINOPT 
LINOPT 
LINOPT 
LINOPT 
LINOPT 
LINOPT 
LINOPT 
LINOPT 
LINOPT 
LINOPT 
LINQF  T 
LINOPT 
LINOPT 
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M  VARIABLES  (SECONDARY  VARIABLES'  IN  TERMS 
OF  THE  REMAINING  N  VAR TABLES  \ PR I MARY 
VARIABLES ) .  K,  A  PESr'JTAII ON  0”  (  1  , . , . , N+rt) , 

SPECIFIES  SUCH  A  PARTITION  INTO  PRIMARY  AND 
5ECGNDAB i  VARIABLES.  K(  I ) , . . .  ,N<M)  ARE  THE 
INDICES  OF  THE  PRIMARY  VARIABLES.  K<M+I 
K(N+N>  ARE  THE  INDICES  OF  THE  SECONDARY 
VARIABLES.  FOR  THE  DUAL  VARIABLES  U'j), 

J  =  1,...,N+H,  PRIMARY  AND  SECONDARY  INDICES 
SWITCH  ROLES,  U(K(N*1 >),... ,UIK(N*H> >  BEING 
PRIMARY .  A  BASIC  SOLUTION  IS  SPECIFIED  BY 
SETTING  EACH  PRIMAL  PRIMARY  VARIABLE  TO 
EITHER  OF  ITS  BOUNDS  AND  EACH  DUAL  PRIMARY 
VARIABLE  TO  ZERO.  THE  INPUT  VALUES  QP  THE 
PRIMAL  PRIMARY  VARIABLES  ARE  SWITCHED  TO  THE 
OPPOSITE  BOUND  IF  NECESSARY  IN  ORDER  TO  CREATE 
A  DUAL-FEASIBLE  BASIC  SOLUTION. 

X  PRIMAL  SOLUTION  ARRAY.  (REAL) 

X(K( J) )  MUST  BE  SET  TO  EITHER  BL(K( J) )  OR 
BU(KU)),  J  *  THESE  ARE  DEFAULT 

VALUES  TO  BE  USED  WHEN  A  VANISHING  U<K< J> ) 

HAKES  X ( K ( J ) >  INDETERMINATE  IN  SETTING  UP 
A  DUAL-FEASIBLE  SOLUTION. 

OUTPUTS— THE  FOLLOWING  VARIABLES  AND  ARRAYS  ARE  DEFINED 
-  OR  REDEFINED  ON  EXIT 


NUMBER  OF  ITERATIONS  SINCE  THE  LAST 
INVERSION. 


(INTEGER) 


ERROR  FLAG.  (INTEGER) 

I  ERR  3  0— OPTIMUM  FOUND. 

1  — INCONSISTENT  CONSTRAINTS. 

2  — ITERATION  LIMIT  REACHED. 

3— INVERSION  FAILED  (BAD  INITIAL  BASIS). 


BASIC  SOLUTION  KEY, 

SET  FOR  THE  CURRENT  BASIS. 

PRIMAL  SOLUTION  ARRAY. 


DUAL  SOLUTION  ARRAY. 

U(J)  IS  THE  DUAL  VARIABLE  (LAGRANGE 
MULTIPLIER)  FOR  THE  CONSTRAINTS  ON  X(J>. 

IT  IS  POSITIVE  IF  THE  UPPER  BOUND  IS  ACTIVE, 
NEGATIVE  IF  THE  LOWER  BOUND  IS  ACTIVE. 

INVERSE  MATRIX  ARRAY. 

X(I)  5  SUM  (J  3  N)  E ( I , J )  *  X(K(J)>, 

I  ~  1  ,  . .  . ,  N . 


(INTEGER) 


(REAL ) 


(REAL) 


(REAL ) 


LIN OPT 
LI, -OPT 
LINQPT 
LliiDPT 
LIN OPT 
LINOPT 
UNOPT 
UNOPT 
UNOPT 
LINQPT 
UNOPT 
IINOPT 
LINQPT 
LINOPT 
UNOPT 
LINOPT 
LINOPT 
LINOPT 
LINOPT 
LINOPT 
LINOPT 
LINQPT 
IINOPT 
LINOPT 
LINOPT 
LINOPT 
LINQPT 
LINOPT 
LINOPT 
LINQPT 
UN  OPT 
LINOPT 
LINOPT 
LINOPT 
UNOPT 
LINQPT 
LINOPT 
LINOPT 
LINOPT 
LINOPT 
LINOPT 
LINOPT 
LINOPT 
LINOPT 
UNOPT 
LINOPT 
LINOPT 
LINOPT 
LINQPT 
LINOPT 
L INOP  ' 
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c 

L I  N  G  P  T 

*  -- 

c 

SCR 

(ALIAS  ARQU  IN  SUBROUTINES)  SCRATCH  ARRA  ( . 

LIHOPT 

1  4  3 

c 

LINGPT 

M6 

c 

LING? T 

1  47 

c 

LINOPT 

143 

c 

NININUrt  DECLARED  ARRAY  SIZES— 

LINO?  T 

1  4? 

c 

LIHOPT 

150 

£ 

A 

ri  f  N 

LINO?  7 

13! 

c 

RDU 

N 

LINOPT 

152 

c 

COL 

N 

LINGPT 

1 53 

c 

31 

H  +  N 

LINOPT 

15-» 

c 

BU 

N  ♦  N 

LINGPT 

155 

c 

K 

H  ♦  N 

LIHOPT 

1  JO 

c 

X 

N  +  N 

LIHOPT 

15? 

c 

U 

N  +  N 

LIHOPT 

153 

c 

E 

H  *  N 

LINOPT 

15? 

c 

SCR 

N 

LII4QPT 

1 40 

c 

LINOPT 

161 

c 

LIHOPT 

162 

c 

LINOPT 

163 

c 

SUBROUTINE 

TABLO  (Q.V.)  PRINTS  THE  FULL  EXPLICIT  TABLEAU. 

LINOPT 

1  o4 

c 

IT  IS  NOT 

CALLED  THROUGH  UNOPT  AND  NU3T  BE  CALLED  SEPARATELY. 

LINOPT 

1  s5 

C  LINE  2 

C - * - - — - — - -  LINE  3 


BINErtSIGN  81(1),  BU(I),  K ( 1 ) ,  X(l),  U(l),  SCRH)  LIHOPT 

C  LINOPT 

COHHGII  /XXXI?/  IOBJ,  H,  N,  HI  N,  INVERT,  ITnrtX,  EPS,  ITER,  I  Erf  R  /XXXI?/ 

COHMGN  /XXXLP /  UP  1,  NPH,  IPIV,  JPIV,  NEGV  /XXXLP/ 

CDH.ION  /XXXLP/  BIGH  /XXXLP/ 

LOGICAL  NIN,  INVERT,  NEGV  /XXXLP/ 

C  LINE 

C . - . - . - . . LINE 

C  LINE 

C  THE  VARIABLE  SIGN  REPRESENTS  A  VERY  LARGE  NL*3ER.  THE  DEFAULT  LINOF'T 

C  VALUE  IS  1  .El  00.  THE  USER  NAY  RESET  THIS  VALLE  IF  SO  DESIRED.  LINCPT 

C  SIGN  OR  -BIGH  NAY  BE  USED  TO  FILL  IN  HISSING  UPPER  OR  LOVER  LINO? 7 

C  BOUNDS.  LINOF'T 

C  LI  NOP T 

DATA  BIGH  /I .El 00/  LIHOPT 

C  LINE 

C - LINE 

C  LINE 

HP  1  =  N  +  t  LIHOPT 

NPH  =  N  +  H  LINGPT 

IF  (.NOT.  INVERT)  GO  TO  10  LIHOPT 

CALL  SE7INV(A,ROU,COl,K,E,5CR>  LIHOPT 

IF  (IERR.E3.3)  RETURN  LINO? T 

10  CONTINUE  .  LINOPT 

CALL  GETROU A,RGU,C0L,E,I08J,SCR)  LINCPT 

DO  50  J  -  ) ,  N  LIhC* I 

H  j  3  K ( J •  L I  HOP' 


16? 

163 

-i 

3 

A 

o 

2 

3 

4 

1?! 

172 

173 

174 
i  73 
1  74 


30 
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If  (iilil)  3CR(J)  =  -SCR < J ) 

Ll.iGPT 

137 

U< KJ)  =  SCR(J) 

LIilOP  f 

las 

If  (U(KJ))  20,  40,  30 

LINOPT 

*  n 

J  'j  r 

NEGATIVE 

LI. '(OPT 

no 

20 

CONTINUE 

LINOPT 

i?i 

X  ( K  J )  =  8KKJ) 

LINOPT 

192 

GO  TO  40 

LINOPT 

193 

POSITIVE 

i.  IN  OPT 

1  ?4 

30 

CONTINUE 

LINOPT 

195 

X(KJ)  =  BU(KJ) 

UNOPT 

196 

40 

CONTINUE 

LINOPT 

197 

50 

CONTINUE 

LINOPT 

193 

DO  60  J  =  NPt ,  NPH 

LINOPT 

199 

U  v  K ( J ) )  =  0. 

LINOPT 

200 

40 

CONTINUE 

LINOPT 

20! 

CALL  PSOL(A,ROU,COL,K,X,E> 

LINOPT 

202 

CALL  DSI«P(A,RQW,COL,BL,BU,K,X,U,E,SCR> 

LINOPT 

203 

If  UERR.NE.O)  GO  TO  110 

LINOPT 

204 

ROUND  X-VALUES  WITHIN  EPS  OF  BOUNDS 

LINOPT 

205 

DO  100  I  =  t,  NPH 

LINOPT 

206 

If  <A6S<X(I)-BL<I)>.LE.£PS)  X(I)  =  BL<I) 

LINOPT 

207 

IF  (A8SUm-BUmJ.LE.EPS)  X(I)  =  BU(I) 

LINOPT 

208 

100 

CONTINUE 

LINOPT 

209 

no 

CONTINUE 

LINOPT 

210 

RETURN 

LINOPT 

21 ; 

END 

LINOPT 

212 

31 
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SUBROUTINE  DSI.1P(A,RQU,C0L,3L,BU 


,  K ,  X  ,!J ,  E  ,  AROU ) 


DUAL  SIMPLEX  ALGORITHM 


DIMENSION  BL(1>,  BU(I),  K<1),  X(l),  U(l),  A K 0 U ( ! ) 

COMMON  /XXX LP/  IOBJ,  M,  N,  MIN,  INVERT,  ITNAX,  EPS,  ITER,  ItRft 
COMMON  /XXXLP/  NP1 ,  NPM,  IPIV,  JPIV,  NEGV 
COMMON  /XXXLP/  8IGM 
LOGICAL  MIN,  INVERT,  NEGV 


IERR  =  2 

UHEN  ITMAX.LT. t  THE  LOOP  IS  PARTIALLY  EXECUTED  '  - 

DO  100  II  =  1,  ITHAX 
CALL  PIVRQU(BL,BU,K,X) 

IF  (IPIV.GT.O)  GO  TO  10 

NO  PIVOT  RQU  INDICATES  THAT  X  IS  OPTIMAL 
IERR  =  0 
RETURN 
10  CONTINUE 

KROU  =«  X(IPIV) 

CALL  GETROWIA, SOU, COL, E, KROU. AROU) 

CALL  PIVCOL(BL,BU,K,X,U,AROU) 

IF  ( JPIV.GT.O)  GO  TO  40 

NO  PIVOT  COLUMN  INDICATES  THAT  THE  CONSTRAINTS 
ARE  INCONSISTENT 
IERR  =  1 
RETURN 
40  CONTINUE 

IF  ( ITMAX.LT . 1 )  RETURN 
C  NEW  SOLUTION  KEY 

KUPIV)  =  K(JPIV) 

K(JPIV)  =  KROU 
CALL  NtUlM(E,AROU> 

C  NEW  DUAL  SOLUTION 

CALL  GcTROU(A,ROU,COL,E,IOBJ, AROU) 

50  70  J  =  I,  H 

IF  (MIN)  AROW(J)  =  -AROW(J) 

U(K( J) )  =  AROU(J) 

70  CONTINUE 

U(K(IPIV»)  =  0. 

C  NEW  PRIMAL  SOLUTION 

X'KROU)  =  BU(KROU) 

IF  ( NEGV)  X(KROU)  =  BL(K'ROU) 

CALL  PSOL > A,ROU,CCL,K,X,E) 

ITER  =  ITER  ♦  I 
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SUBROUTINE  ?1VRQUUI,BU,K,X) 
C 

C . 

PIVOT  ROW  SELECTION 


DIMENSION  Bid),  BUM),  KU),  XU) 

COMMON  /XXX IP/  ItiBJ,  M,  N,  MIN,  INVERT,  ITHAX,  EPS,  ITER,  ItRR 
COMMON  am.?/  NPf ,  NPN,  IPIV,  JPIV,  NEGV 
COMMON  /XXXLP/  BIGM 
LOGICAL  MIN,  INVERT,  NEGV 


I?IV  *  0 

If  (NPH.LT.NP1)  RETURN 
VIOL  =  0. 

DO  50  II  =  NP1 ,  NPH 
I  =  Kill) 

C  CHECK  CONSTRAINTS  ON  XU) 

D  =  XII)  -  BLU) 

IF  (D.GE.-EFS)  GO  TO  10 
0  =  -B 

IF  (VIOl.GT.D)  GO  TO  40 
VIOL  *  D 
IPIV  =  II 
NEGV  *  .TRUE. 

GO  TO  40 

10  CONTINUE 

B  =  XU)  -  sum 
IF  (D.LE.EPS)  GO  TO  30 
IF  (VIOL.GT.D)  GO  TO  40 
VIOL  =  B 
IPIV  =  II 
NEGV  FALSE. 

30  CONTINUE 

40  CONTINUE 

30  CONTINUE 

RETURN 
£NB 
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SUBROUTINE  PIVCQL(BL,B«,K,Xt«.flROtf)  PIVCOL 

LI, HE 

- 

LI  WE 

PIVOT  COLUMN  SELECTION  PIVCOL 

LINE 

— . - . - . - . —  LINE 

LINE 

DIMENSION  Bl(l),  3U<1>,  KU ) ,  X(1),  UC 1 ) ,  AROU(I)  PIVCOL 

PIVCOL 

COMMON  /XXX LP/  IOBJ,  H,  N,  MIN,  INVERT,  UMAX,  EPS,  ITER,  I ERR  /XXXLP/ 

COMMON  /XXXLP/  N PI,  NPrt.  IPIV,  JPIV,  NEOV  /XXXLP/ 

COMMON  /XXXLP/  BIGM  /XXXLP/ 

LOGICAL  MIN,  INVERT,  NEGV  /XXXLP/ 

LINE 

. . LINE 

LINE 

JPIV  =  0  PIVCOL 

ii  =  BIGM  PIVCOL 

DO  30  JJ  =  1,  N  PIVCOL 

J  =  K(JJ)  PIVCOL 

AA  =  AROU(JJ)  PIVCOL 

IF  (NEGV)  AA  -  -AA  PIVCOL 

IF  < AA.GE.O.  .AND.  X( J) .EQ.BK J) )  GO  TO  20  PIVCOL 

IF  (AA.LE.0.  .AND.  X C J ) .EQ.BU< J) )  GO  TO  20  PIVCOL 

R  =  U(J)/AA  PIVCOL 

IF  (R.GT.U)  GO  TO  10  PIVCOL 

IF  (R.EQ.U  .AND.  RANF(AA) .uT.0.5)  GO  TO  10  PIVCOL 

JPIV  =  JJ  PIVCOL 

U  =  R  PIVCOL 

10  CONTINUE  PIVCOL 

20  CONTINUE  PIVCOL 

30  CONTINUE  PIVCOL 

RETURN  PIVCOL 

END  PIVCOL 
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C 


SUBROUTINE  GETR0U(A,R0U,C0L,E,KR0U,ARGU) 


GENERATION  OF  CONSTRAINT  COEFFICIENTS  FOR  THE  CURRENT  BASIS 


DIMENSION  A(l),  ROUd),  CQld),  Ed),  AROWd) 

INTEGER  RQU,  COL 

COMMON  /XXX LP/  IOBJ,  M,  N,  MIN,  INVERT,  UMAX,  EPS,  ITER,  TERR 
COMMON  /XXXLP/  NP1 ,  NF'M,  IPIV,  JPIV,  NEGV 
COMMON  /XXXLP/  BI6H 
LOGICAL  MIN,  INVERT,  NEGV 


IF  (XRQU.GT.N)  GO  TO  20 

ORIGINAL  INDEPENDENT  VARIABLE.  GET  RON  KRO’4  OF  THE  INVERSE. 
JJ  =  0 

DC  10  J  *  1,  N 

AROU(J)  =  E(KROU+ JJ) 

IF  (ABS< AROU( J) ) . LE . EPS )  ARQU(J)  =  0. 

JJ  =  JJ  +  N 
10  CONTINUE 

GO  TO  50 
20  CONTINUE 

ORIGINAL  DEPENDENT  VARIABLE. 

MULTIPLY  ORIGINAL  RQU  BY  THE  INVERSE. 

KK  =  ROU(KROU-M) 

JJ  *  0 

DO  40  J  *  t,  N 
AROU(J)  =  0. 

DO  30  I  =  1,  N 

AROUiJ)  =  AROU(J)  ♦  A(KK+COL(I))*E(I+JJ) 

30  CONTINUE 

IF  (ABS( ARUU( J) ) .LE.EPS)  AROU(J)  =  0. 

JJ  =  JJ  +  N 
40  CONTINUE 

50  CONTINUE 
RETURN 
END 
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SUBROUTINE  FSQL(A,RQV,CGL,X,X,E) 
C 

C— . 

c 

C  PRIMAL  SOLUTION 


DIMENSION  A(1),  ROUd),  COLd),  Kd),  XI 1  > ,  E(I) 

INTEGER  ROW,  COL 

COMMON  /XXX LP/  IOBJ,  M,  N,  HIN,  INVERT,  UMAX,  EPS,  ITER,  I  ERR 
COMMON  /XXX LP/  NPI ,  NPM,  IPIV,  JPIV,  NEGV 
COMMON  /XXX LP/  BIGM 
LOGICAL  MIN,  INVERT,  NEGV 


DO  30  I  =  NPI ,  NPM 
KI  =  K( I ) 

IF  (KI.GT.N)  GO  TO  20 
X(KI)  =  0. 

JJ  =  0 

DO  10  J  =  1,  N 

X ( K I >  =  X ( KI )  +  EIKI+JJ)  +  X ( K  <  J ) ) 

JJ  =  JJ  +  N 
10  CONTINUE 

20  CONTINUE 

30  CONTINUE 

DO  60  I  =  NPI ,  NPM 
KI  =  KID 

IF  (KI.LE.N)  GO  TO  50 
X ( K I )  =  0. 

KK  =  ROU(KI-N) 

00  40  J  =  1,  N 

X ( K I )  =  X (KI >  +  A(KK+COL( J) )  *  X ( J) 
40  CONTINUE 

50  CONTINUE 

60  CONTINUE 

RETURN 
END 
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SUBROUTINE  SETINV!  A,  ROU,  COL,  K,Z,ARQ!J} 


INITIAL  INVERSE 


DIMENSION  K  ( I ) ,  Em,  ARQU(I) 

COMMON  / X X XLP /  IOBJ,  M,  N,  MIN,  INVERT,  UMAX,  EPS,  ITER,  IERR 
COMMON  /XXXLP/  NP1 ,  NPH,  IPIV,  JPIV,  NEGV 
COMMON  /XXXLP/  BIuH 
LOGICAL  MIN,  INVERT,  NEGV 


SET  E  TO  THE  IDENTITY 
JJ  =  0 

DO  20  J  =  1,  N 
DO  10  I  -  !,  N 
EU+JJ)  =  0. 

10  CONTINUE 

EU+JJ)  =  1. 

JJ  =  JJ  ♦  « 

20  CONTINUE 

C  GENERATE  INITIAL  INVERSE 
DO  30  J  =  1,  N 
K(J)  =  — K  <  J ) 

30  CONTINUE 

DO  90  JJ  -  1 ,  N 

DO  40  J  =  1,  N 

IF  (K(J).LT.O)  GO  TO  50 

40  CONTINUE 

50  CONTINUE 

KROU  =  -K( J) 

CALL  GETROU(A,ROU,CGL,E,KRQU, AROU) 

ROUMAX  =  0. 

DO  70  L  =  1,  N 

TEST  =  ABS( AROU(L) ) 

IF  (K(L).GT.O  .OR.  TEST. LT. ROUMAX)  GO  TO  30 
ROUMAX  =  TEST 
JPIV  =  L 

60  CONTINUE 

'■70  CONTINUE 

IF  (ROUMAX. GT.O.)  GO  TO  30 

IERR  =  3 
RETURN 
30  CONTINUE 
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90 


K(J)  =  K(JPIV) 

3£7Ii'<‘> 

42 

K(JPIV)  =  KROU 

SET INU 

43 

CALL  NEUINV(E,ARQU) 

SET I  MV 

44 

CONTINUE 

St  UNO 
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ITER  =  0 

StTINV 
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RETURN 
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c 

n- 

SUBROUTINE  NEUINVIE, AROU) 

NZJINV 

LINE 

n 

n 

U 

Ll  nil 

5 

c 

LINE 

4 

c 

INVERSE  UPDATE  BY  COLUMN  OPERATIONS 

NEUINV 

4 

c 

LINE 

2 

L  1  >  T  C 

3 

c 

LINE 

4 

DIMENSION  Ed),  AROUd) 

NEUINV 

6 

c 

NEUINV 

7 

COMMON  /XXXLP/  IOBJ,  M,  N,  MIN,  INVERT,  UMAX,  EPS,  ITER,  I ERR 

/XXXLP/ 

2 

COMMON  /XXXLP/  NP),  NPN,  IPIV,  JPIV,  NEGV 

/XXXLP/ 

3 

COMMON  /XXXLP/  BIGM 

/XXXLP/ 

4 

LGGICAL  MIN,  INVERT,  NEGV 

/XXXLP/ 

5 

c 

LINE 

2 

r 

3 

U 

—  lltit 

c 

LINE 

4 

JJPIV  =  <  JPIV— 1 )*N 

NEUINV 

10 

DO  20  I  »  1,  N 

NEUINV 

11 

EPIV  -  E<I+JJPIV)/AROU(JPIV) 

NEUINV 

12 

JJ  =  0 

NEUINV 

1 3 

DO  10  J  =  I,  N 

NEUINV 

14 

E( I+JJ)  =  Eil+JJ)  -  EPIV*AROU( J) 

NEUINV 

15 

JJ  =  JJ  +  N 

NEUINV 

16 

10 

CONTINUE 

NEUINV 

17 

E( I+JJPIV)  =  EPIV 

NEUINV 

IS 

20 

CONTINUE 

NEUINV 

19 

RETURN 

NEUINV 

20 

END 

NEUINV 
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SUBROUTINE  TA?LO<A,ROU,COL, BL,BU,K. X,E, SCR, KORD; 


TABLEAU  PRINTOUT 

KORD  IS  AN  ARRAY  OF  LENGTH  AT  LEAST  N  USED  FOR  REORDERING 
/(!),...  ,Kin)  IN  ASCENDING  ORDER. 

SEE  LINOPT  FOR  DESCRIPTIONS  OF  THE  OTHER  PARAMETERS. 

LINQPT  MUST  HAVE  BEEN  CALLED  BEFORE  CALLING  TABLO. 


DIMENSION  AM  ) ,  ROUM),  COLM),  SL(!>,  BUM  >,  KM),  XM),  EM) 
DIMENSION  SCR  \  1  > ,  KORDM) 

INTEGER  RGU,  COL 

COMMON  /XXXLP/  IOBJ,  M,  N,  MIN,  INVERT,  ITMAX,  EPS,  ITER,  IERR 
COMMON  /XXXLP/  NP1 ,  NPM,  IPIV,  JPIV,  NEGV 
COMMON  /XXXLP/  BIGM 
LOGICAL  MIN.  INVERT,  NEGV 


1  FORMAT  MHI//T55,*T  A  B  L  E  A  U*//tX,*  I  *,5X,*BL(I)*, 

/  6X,*X<I)*,5X,*BU(I)*,1X,10A10/(36X,10A10>) 

2  FORMAT  MHO,I3,)X,3F10.3,1X,1OF1O.3/(36X,1OF1Q.3)) 

3  FORMAT  (2X,*T<I,*,I3,*)*> 

4  FORMAT  (///IHO,*THE  OBJECTIVE  VARIABLE,  X(*,I3,*>,  IS  TO  BE  *,A10) 


DO  MO  J  =  1,  H 
KORD ( J )  =  J 
MO  CONTINUE 
DO  130  J  =  1,  N 
JNIN  =  J 

DO  120  JJ  =  J,  N 

IF  (K(KORD(JJ)I.LT.K(KORD(JMIN)))  JMIN  =  JJ 
120  CONTINUE 

KTEiiF  =  KORD'.  J) 
kORD(J)  =  KORD(JMIN) 

KORDTJMIN)  =  KTEMP 
130  CONTINUE 

DO  1 0  J  =  1 ,  N 

ENCODE  (10,3,  SCR  <  J) )  K(KORDU)) 

10  CONTINUE 

PRINT  1,  (SCR(J),  J  =  1,  N) 
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DO  20  I  =  1 ,  NPii 

TABLO 

40 

CALI  GETR0U(A,R0U,C0L,E,I,SCR) 

T  ABLD 

41 

PRINT  2,  I.  3L(I),  X(I),  BU(I),  ( SCR ( KORD ( j > ) ,  J  =  1,  N) 

TABLO 

42 

20  CONTINUE 

TABLO 

43 

OPT  =  IOHWAXIMIZED. 

TABLO 

44 

IF  (HIM)  OPT  =  1 0HM I H I M I ZED . 

TABLO 

45 

PRINT  4,  IC3J ,  OPT 

7A8L0 

4.j 

RETURN 

TA3LQ 

47 

END 

TABLO 

43 

42 
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